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A non-linear conjugate gradient optimization scheme is used to obtain excitation energies within 
the Random Phase Approximation (RPA). The solutions to the RPA eigenvalue equation are located 
through a variational characterization using a modified Thouless functional, which is based upon 
an asymmetric Rayleigh quotient, in an orthogonalized atomic orbital representation. In this way, 
the computational bottleneck of calculating molecular orbitals is avoided. The variational space 
is reduced to the physically-relevant transitions by projections. The feasibility of an RPA imple- 
mentation scaling linearly with system size, N, is investigated by monitoring convergence behavior 
with respect to the quality of initial guess and sensitivity to noise under thresholding, both for 
well- and ill-conditioned problems. The molecular-orbital-free algorithm is found to be robust and 
computationally efficient providing a first step toward a large-scale, reduced complexity calcula- 
tion of time-dependent optical properties and linear response. The algorithm is extensible to other 
forms of time-dependent perturbation theory including, but not limited to, time-dependent Density 
Functional theory. 
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I. INTRODUCTION 

Matter responds to electromagnetic perturbation in a 
time- dependent fashion: incident light induces periodic 
fluctuations within the electron density of a molecule that 
can be described by its excitation spectrum. The exci- 
tation spectrum is of fundamental importance to many 
fields, ranging from analysis of interstellar clouds to the 
molecular basis of disease. Unfortunately, excited states 
arc difficult to calculate for large complex systems be- 
cause the scaling of the computational cost with the num- 
ber of atoms is prohibitive. While numerous efforts have 
been devoted to the development of reduced complex- 
ity algorithms for ground state properties, 1 much less 
work has been focused on efficient algorithms for excited 
state response properties. The purpose of this paper is 
to investigate a method for variational characterization 
of the excitation spectrum that could potentially scale 
linearly with system size. This would allow studies of 
much larger systems than currently achievable. The ex- 
citation spectrum is described by the Random Phase Ap- 
proximation (RPA) within time-dependent Hartree-Fock 
theory, 2,3 ' 4 ' 5 ' 6 ' 7,8 but our algorithm is general and can 
be applied also to time-dependent Density Functional 
Theory. 9 ' 10 



A. The RPA Equation 

Concomitant to the development of Many Body the- 
ory to describe the ground states of molecules, work to 
calculate properties of the more elusive excited states em- 
ploying the Random Phase Approximation (RPA) began 
in the early 1950's. Avoiding the complications of ad- 
dressing independent particles in a many electron sys- 



tem, the original, classical mechanical RPA treats the 
electron-repulsion terms as part of an ensemble average. 
The Fourier transforms of the Coulomb terms have "ran- 
dom phases" that cancel, hence the name. 11,12 Recog- 
nition of an explicit quantum mechanical connection to 
single determinants 13 led to the demonstration of equiv- 
alence between the RPA and a time-dependent extension 
of Hartree-Fock (HF) theory - thus permitting a fully 
quantum mechanical treatment of matter under light- 
induced perturbation. 14 

During the 1960's, three equivalent formalisms de- 
veloped around application of the RPA to calculate 
excited states. While these derivations, based upon 
Equations of Motion, 15,16 ' 17 Green's functions, 18,19 and 
time-dependent Hartree Fock theory 20 ' 21,22 are non- 
trivial, 6,7,23 it is sufficient to note that for electronic tran- 
sitions, a description utilizing only the one-body den- 
sity is valid - provided that the particle excitation en- 
ergies arc smaller than the Fermi energy and two-body 
correlations can be neglected. 20 In this case, electronic 
excitations are well-described by the RPA eigenvalue 
equation, 6,23 written in the Molecular Orbital, "MO" ba- 
sis as 

The resonant frequencies, or excitation energies, are rep- 
resented by the eigenvalues, u>. The elements of the ma- 
trices A and B are given by 

Ami,nj — (^m €i)&ij&mn ~t~ Vmj.in Vmj,ni : (2) 

and 

Bmi,nj — Vmn,ij Vmnji > (3) 
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and the elements of V mn ,ij are the conventional two- 
electron integrals. 8 The i, j indices are from the set of oc- 
cupied states while the m, n indices correspond to the vir- 
tual orbitals, while e m and e, denote the Fockian eigenen- 
ergies. 

The matrices A and B correspond to 4th order tensors 
of dimension (iV occ x N vlrt ) x (N virt x N occ ) spanning 
the Liouville space of transitions between the occupied 
(occ) and virtual (virt) subspaces. These act upon the 
vectors X and Y, composed of orbital coefficients, so that 
particle-hole (ph) transitions arc described by X while Y 
contains the hole-particle (hp) transitions. 

The first term of A in Eq. (2) corresponds to the un- 
dressed, bare excitations, i.e., those predicted by Koop- 
man's theorem. The last two terms, or B, add a 
correlation-based correction to the bare energies of A 
based upon the Coulomb and exchange interactions. 
(Setting Y = 0, produces the Tamm-Dancoff 24, 25,26 ap- 
proximation.) Finally, the unitary matrix, IN = (q _?/), 
is a unit diagonal metric tensor, serving as an orthonor- 
malization constraint, 19,27 defining the indefinite inner 
product associated with the space of ph-hp transitions, 
in the Molecular Orbital basis. 



B. Linear Scaling Approaches to Solving the RPA 
Equation 

The RPA equation was originally derived in the molec- 
ular orbital representation, as in Eq. (1), and the fa- 
miliarity of "molecular orbitals" in discussions involv- 
ing ground states render it a popular basis in which 
to work. 10 ' 28 However, the molecular-orbital representa- 
tion requires a full eigenfunction solution of the ground 
state problem, which typically requires a computational 
cost that scales as 0(N 3 ), where N is the number of 
basis-functions, assumed to be proportional to system 
size. A requirement for any reduction of this computa- 
tional 0(N 3 ) bottleneck is, therefore, to find a molecular- 
orbital-free algorithm for the solution of the RPA equa- 
tion. 

Recently, a number of groups 29 ' 30,31,32,33,34 have 
achieved a linear scaling computational complexity for 
the ground state self-consistent field (SCF) problem in 
Hartree-Fock (or density functional theory) using "fast" 
algorithms for computation of the Fockian F and sparse 
matrix algebra (dropping of small elements) to exploit 
quantum locality of the density matrix P. If the tran- 
sition densities in the time-dependent response equa- 
tions also demonstrates quantum locality, then the same 
fast methods used for the ground state problem are 
applicable. 30,33 

Solving the time-dependent quantum response prob- 
lem in O(N) is pivotal in studies of large scale systems 
currently inaccessible to conventional methods. Perhaps 
the most successful approach, to date, is to propagate an 
electron impulse response through numerical integration 
in real time, 35,36,37,38 and then retrieve the spectra from 



the time series through Fourier transformation. More 
recently, Coriani, et al. 29 have implemented the matrix 
exponential approach of Larsen, et al., 34 observing an ac- 
celeration in computation of excitation energies for one 
dimensional systems. 

Another reduced complexity approach for time- 
dependent response calculations was recently presented 
by Izmaylov, et al. 30 and Kussman and Ochscnficld. 39 
It is worth noting that in the adiabatic zero-frequency 
limit, when uj — > 0, the adiabatic response problem can 
be solved with surprising efficiency in linear scaling com- 
plexity using adiabatic density matrix perturbation the- 
ory based on purification. 40 Linear scaling density matrix 
perturbation theory can be applied to the calculation of 
response properties of molecules, both for lower 41 and 
higher order perturbations, 32 as well as for the crystalline 
problem, 31 including the electric polarizability. 

The reduced complexity approach in this paper is 
based upon a well-established variational characteriza- 
tion of the eigenvalue spectrum as applied to the RPA 
equation. The key idea is to use an molecular-orbital- 
free approach, avoiding the 0(N 3 ), bottleneck. This 
is achievable through a functional optimization of an 
asymmetric Raylcigh Quotient as formulated by Thou- 
less more than four decades ago. 19 The intent of this 
paper is not to present a linear scaling algorithm, but 
to analyze and discuss the limitations and feasibility of a 
variational optimization of a Thouless functional in the 
context of reduced complexity calculations. 

II. MOLECULAR-ORBITAL-FREE 
TIME-DEPENDENT PERTURBATION THEORY 

To derive a molecular-orbital-free formulation for the 
RPA equation suitable to O(N) calculations, we may 
start from time-dependent Hartree-Fock theory, 6,23 

= [F, P] s = FPS - SPF, (4) 

where S is the overlap matrix, P is the single-particle 
density matrix of the Hartree-Fock ground state, and F 
is the effective single-particle Hamiltonian, i.e., the Fock- 
ian (or the Kohn-Sham Hamilitonian, in a generalization 
to time-dependent Density Functional theory) . In an or- 
thogonalizcd representation, S becomes the identity ma- 
trix. 

Looking at the first-order response under variation of 
the density matrix SP, we find that, 

dSP 

i— = [F,6P] s + [G(6P),P} s , (5) 

which, in the frequency domain, gives the RPA linear 
response eigenvalue equation, 

[F,x] B + [G(x),P] a =ux. (6) 

Here, x is the Fourier transform of SP and 

G(x)=2J[x)-K[x}. (7) 
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The left commutator in Eq. (6) gives the zeroth-order 
approximation corresponding the bare excitations, and 
the second commutator with G(x), includes additional 
Coulomb, J[x], and exchange screening, If [a;]. In a gen- 
eralization to time-dependent DFT, the exchange screen- 
ing is replaced by the exchange correlation screening, i.e., 
the second-order functional derivative of the exchange- 
correlation action. 9,42 ' 43,44 The RPA excitation spectrum 
is thus given by the cigcnfrequcncics lo corresponding to 
ph,hp transitions in Eq. (6). 

In a compact form, we can express the RPA equation, 
Eq. (6) as 



Lx 



LOX. 



(8) 



The vector x is dyadic, corresponding to the unrolled 
N x N matrix a?, i.e., x NxN -w* Xjy2 xl , where the double- 
headed arrow denotes both equivalence and a tensorial 
mapping, or simply a stack operation. 45 For the matrix 
transpose x T , we use the corresponding unrolled vector 
notation x l . In the following, we employ a mixed super- 
vector/matrix notation; projection is most natural for 
matrices, while the use of a vector notation lends itself 
to gradient-based minimization. The action of E onto x 
in Eq.(8) is thus given by 



Ex & [F,x] s + [G(x),P] 



(9) 



The general formulation of the RPA equation, as ex- 
pressed in Eqs. (6) and (8), is independent of the basis- 
set representation and can thus be applied in a molccular- 
orbital-frcc approach, avoiding an expensive diagonal- 
ization scaling as 0(N 3 ) with system size N. In our 
molecular-orbital-frcc algorithm, we employ an orthogo- 
nalizcd atomic orbital basis representation, with S = I, 
which can be efficiently constructed with O(N) complex- 
ity through a congruence transformation, 46 e.g., based 
upon the approximate inverse Cholesky transform. 47 
While this work involves solution of the RPA eigenvalue 
equation, extension to time-dependent Density Func- 
tional theory is straightforward. 



III. NON-LINEAR CONJUGATE GRADIENT 
OPTIMIZATION OF THE THOULESS 
FUNCTIONAL 

While the Lanczos algorithm has been used to itera- 
tively solve the RPA equation, 48 severe problems are ex- 
perienced when calculating high- lying excitations, e.g., 
due to orthogonality constraints. 49,50 More importantly, 
achieving linear scaling complexity requires sparse lin- 
ear algebra, which may preclude the Lanczos algorithm 
due to numerical instabilities. 51,52 ' 53,54 Our molecular- 
orbital-free scheme utilizes a non-linear conjugate gradi- 
ent optimization of a Raylcigh quotient related to the 
method of Muta. 55 The use of non-linear conjugate gra- 
dients are particularly advantageous in the context of 



linear scaling algorithms, because of the ability to re- 
main robust under an incomplete sparse matrix algebra, 
as demonstrated in the work of Simoncini 52 and Notay. 56 
The core of our algorithm is a variational characteri- 
zation of the excitation spectrum based on the Thouless 
functional. 19 Thouless demonstrated the possibility of a 
variational approach to solving the RPA equation via 
iterative optimization of an asymmetric Raylcigh Quo- 
tient. This functional, when expressed in representation- 
independent form, becomes 



n[x] 



x* ■ lux* 



(10) 



where it is understood that the numerator is computed as 
in Eq. (9) and the denominator is given by the absolute 
value of the Euclidean vector product denoted by the 
dot between x and x„. The metric tensor IN is included 
through 



UNTf <s> (P-Q)x 



(11) 



where the subscript IN denotes the action of IN onto x, P 
is the occupied subspace projector and Q = I — P is the 
complimentary projector for the virtual subspace. 

Only stationary solutions to the Thouless functional in 
Eq. (10), corresponding to ph-hp transitions between the 
occupied and virtual subspaces, are of physical relevance. 
Rather than impose ph-hp symmetry explicitly by con- 
struction, it is straightforward to reduce the variational 
search space to the physically-relevant solutions by the 
projection. 7,57 ' 58,59 



or, equivalently, 



x P = PxQ + QxP 



[x,P],P] 



(12) 



(13) 



This projection conforms to the ph-hp formalism of the 
RPA Equation in the MO basis, Eq. (1), with removal of 
non-physical states and reduces the size of the variational 
search space considerably. 

For large, sparse problems, it is possible to con- 
struct the P and Q projectors (or density matrices) 
with linear scaling complexity using recursive purifica- 
tion methods. 60,61,62 In the general case, the metric ten- 
sor IN, which occurs implicitly in the Thouless functional, 
corresponds to the indefinite scalar product 63,64 

(v, u) P & Tx{< K. P}} = Tr{<(P - Q)u p } (14) 

with the norm 



\\x\\ P = \\x\\ P = sJ\{x,x) P \ . 



(15) 



We may now consider optimization of the representation- 
independent Thouless functional, 



Q\x] 



(.r, E.r)i 



(16) 
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This formulation of the Thoulcss functional implicitly in- 
vokes ph-hp symmetry for the excitations, which pro- 
duces paired eigenvalues, ±Wj as would be associated 
with the Liouville operator. 65 

To locate the first transition, (i = 1), we minimize 
Eq. (16), 



Vi = argminf2[x] 

X 

which yields the cigcnfrequcncy 



(17) 



(18) 



The search for subsequent eigenstates requires that lower 
lying eigenvectors be either projected out or shifted away, 
such that they are not rediscovered by consecutive min- 
imization. We use a Wilkinson shift 66 of the interior 
eigenvalues, shifting toj to u>j + a, which is outside the 
region of interest, written as the shifted ~Lx: 



L x + (ujj + a) [<0j (Vj , x) P + Vj (Vj , x)f 

3 



(19) 



where {vj} and {wj} are previously determined RPA 
eigenstates and excitation energies, respectively. 

Our molecular-orbital-free algorithm utilizes a con- 
ventional Polak-Ribierc nonlinear conjugate gradient 
algorithm, 67,68 with restarts, in an orthogonal atomic or- 
bital (e.g., Lowdin 69 ) basis, as is summarized in Fig. 1. 
The outer loop runs over the first M interior eigenvalues, 
while the inner loop iterates over the non-linear conju- 
gate gradient steps. The projections in Lines 4 and 10 
eliminate non-physical states by imposing ph, hp symme- 
try, which significantly reduces the search space. Note: 
the action of the E operator is performed through the 
Fock builds in Eq. (9). The first E operation on Line 
7 is used for calculation of the gradient, and the second 
E operation occurs at Line 13, which is used in the line 
search of Line 14. 

The Wilkinson shift occurs on Lines 6 and 7 while con- 
struction of the Thoulcss functional occurs in Line 8. The 
gradient g, is defined in Lines 9 and 10 and the the con- 
jugate gradient search directions p k are given by the sub- 
sequent calculations on Lines 11 and 12. Restarts for 
(3 = were not necessary, and did not occur during our 
test calculations. 

The functional minimum of the line search along the 
conjugate gradient directions on Line 14 is given by 



A+ = 



-b± s/W- 



Aac 



2ac 



(20) 



where 



a = (p,t) P [ip,x)p + (x,p) P ] 

-{p,p) P [{p,s) P + {x,?) P ] (21) 



FIG. 1: The Molecular-Orbital-Free RPA Algorithm. 

1 for i = 1 to M do 

2 Generate initial vector, x 

3 for k — 1, until convergence 

4 x = VxQ + QxP 



i-l 



6 (uj + a)) [vj(vj,x) P + v-(v-,x)j 



s — Ix + s 
(x,s) P 



n 



x 



9 gk = 2s — 2Q,x 

10 g k =Pg k Q + Qg k P 



j3 = max { — — — * — ' ' ,0 



12 p k = g k + (3p k 

13 t= JLp k 

H Afc = argmin Q[x + Xp k 

x 

15 x — x + Xkf) k 

16 end do 



2{p,t)p{x,x) P - 2(x,s)p(p,p), 



(22) 



17 Vi — X 

18 LUi = ii[Vi] 

19 enddo 



c = (x, X) P [(p, s) P + (x, t) P ] 

- {x,s)p [ip,x) P + (x,p)p] . (23) 

After each inner loop iteration over fc, the desired ith 
cigenpair composed of eigenvector the Hi and eigenvalue 
LUi, are given on Lines 17 and 18. 



IV. PERFORMANCE OF THE 
MOLECULAR-ORBITAL-FREE ALGORITHM 

A. Illustration of the RPA Eigenvalue Spectrum 

To illustrate the performance of the molccular-orbital- 
free solution of the RPA equation, the properties of the 
solutions and various relevant concepts, a schematic pic- 
ture containing a hypothetical set of spectra is provided 
in Fig. 2. 

The spectrum to the extreme left, labeled (FULL), de- 
picts a complete eigenvalue spectrum. All eigenvalues, 
physical and non-physical are included: no projections 
have been performed. There are numerous "bands" that 
might imply clustering or degeneracies which would slow 
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down a variational search of the eigenstates. 

The next spectrum, to the immediate right, (TRAN- 
SITIONS), demonstrates the effect of the projection in 
Eq. (12) including only the subspace of ph,hp transitions. 
The removal of unphysical states significantly reduces the 
density of eigenvalues, particularly around zero. This 
strongly facilitates the search for eigensolutions of the 
RPA equation, in particular for low-lowing excitations of 
physical interest. 70 ' 71,72 ' 73,74 




-600- 



-800 - \ 

FULL TRANSITIONS CALCULATED KOOPMAN'S 

FIG. 2: Illustration of a set of eigenvalue spectra for the RPA 
equation. The entire spectrum, (FULL), is shown on the far- 
left. Notice the dense clustering of eigenvalues around zero. 
Immediately to the right, (TRANSITIONS), depicts the spec- 
trum containing only physical particle-hole and hole-particle 
excitations. The next spectrum, (CALCULATED), is an ex- 
pansion emphasizing a few low-lying transitions. The spec- 
trum to the extreme right, (KOOPMAN's) contains the corre- 
sponding transitions as would be estimated from Koopman's 
theorem. 

The third and fourth spectra, on the right, CALCU- 
LATED and KOOPMANS, provide enlargements of the 
eigenvalue spectrum, in the positive region near zero. 
The CALCULATED spectrum, shows low-lying ph-hp 
transitions, given by the RPA eigenvalue equation. Note 
that high energy excitations are not expected to be well- 
described by the RPA equation. 23,75,76,77 The Koopman's 
spectrum shows the bare excitations, without Coulomb 
or exchange screening, and corresponds to the eigenvalue 
excitation spectrum of the ground state Fockian. Koop- 
man's theorem provides only an approximate spectrum, 8 
but, as will be shown below, initial vectors derived from 
Koopman's theorem produce good initial guesses for our 
non-linear conjugate gradient optimization. 



B. Convergence Behavior 

To study the behavior of our molecular orbital-frce al- 
gorithm, it was prototyped in an orthogonalizcd atomic 



orbital basis employing dense linear algebra and a con- 
ventional 0(N 4 ) approach to Hartree-Fock theory. Our 
test calculations are thus not performed with any lin- 
ear scaling complexity. The present implementation is 
limited to s-type STO-3G basis functions and is not ex- 
pected to produce chemically-relevant data. Rather, we 
utilize this description to characterize the most impor- 
tant features related to more accurate representations, 
as well as consider the problems inherent to linear scal- 
ing implementation. In this manner, we can simulate the 
influence of large basis sets, extended periodic systems 
and complex molecules. 

As with more conventional approaches to solving the 
excitation problem within time-dependent perturbation 
theory, the work required to resolve an eigenstate in- 
creases with an increasing condition number, k, i.e., the 
ratio between the highest and lowest singular value of 
L. To study convergence and other properties, we gen- 
erated test systems using a linear arrangement of four- 
teen hydrogen atoms. Progressively smaller inter-atomic 
spacings produce correspondingly higher condition num- 
bers, k, which parallels that which could arise as the size 
of system or basis set increases. Notice, a variation of 
the condition number also illustrates the effect of pre- 
conditioning. All condition numbers, k, are estimated 
approximately in order of their magnitude. 

To simulate the effects of an incomplete, sparse matrix 
algebra - an absolute necessity for linear scaling capabil- 
ity - we add a random matrix with elements of amplitude 
±t after each application of 1L to a vector, i.e., the La; 
Fock builds. This is equivalent to using a looser numeri- 
cal threshold in the case of a vanishing difference density 
matrix, which has been shown to yield linear scaling. 78 

Two types of initial guesses are considered; a random 
guess, and a "Koopmans' guess." The Koopmans' guess 
was based on direct diagonalization of the ground state 
Fockian. This is an expensive procedure that certainly 
does not scale linearly with system size. The purpose is to 
study the effect of an improved initial guess. An efficient 
0(N) construction of an accurate initial guess remains 
a very important, yet unsolved problem not discussed in 
this article. 79,80,81,82 

The convergence is measured in terms of the relative 
errors of the approximate RPA eigenvalues, e„, where the 
error, 

Err(n) = | £ "~^ / |, (24) 

En 

is calculated in each conjugate gradient iteration, n. The 
reference eigenvalues, 0J ref , were obtained from direct di- 
agonalization of the Ij matrices using the ZGEEV routine 
in the LAPACK library. 83 

1. Small k, Varied Initial Guess 

The convergence behavior for the first 5 eigenvalues, 
Ai — A5, of the system with a condition number k = 10 2 
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is depicted in Fig. 3. The algorithm starts with ran- 
domly generated initial vectors for each of the 5 eigen- 
values sought. The non-linear conjugate gradient opti- 
mization for each eigenvalue is then allowed to proceed 
until Err(rc) < lCT 12 . 



Condition Number, k=10 
Initial Guess: Random Vector 
Eigenvalue Convergence 
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— Eigenvalue 1 

— Eigenvalue 2 
■ — Eigenvalue 3 
Eigenvalue 4 

Eigenvalue 5 



80 100 120 140 160 180 200 
Iteration (n) 

FIG. 3: Convergence of the first 5 eigenvalues for a well- 
conditioned (k = 10 ) matrix using random initial guesses. 



Contrast these results with the curves in Fig. 4, which 
contains the convergence patterns for the same system, 
but this time starting with initial guess vectors based 
upon Koopman's theorem. Despite the apparent non- 
ideality implied by the hypothetical Koopman's-based 
spectrum (KOOPMAN's) in Fig. 2, the initial eigenvec- 
tors generated by Koopman's guess provide significant 
improvement in the rate of convergence compared to that 
of the random vectors used to generate the plots in Fig. 3. 
For the third eigenvalue, the convergence is improved by 
almost an order of magnitude. 



2. Convergence: Varied k, Good Initial Guess 
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Condition Number, k=10 
Initial Guess: Koopman's Theorem 
Eigenvalue Convergence 
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FIG. 4: Convergence of the first 5 eigenvalues for a well- 
conditioned (ft = 10 2 ) matrix using initial guesses derived 
from Koopman's theorem. 
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FIG. 5: Convergence of the first eigenvalue for matrices with 
K € [10 2 , 10 4 ]. In each case the Koopmans' guess was used. 



While the behavior of the algorithm for well- 
conditioned matrices is useful for proof-of-concept, the 
performance of any algorithm in the presence of ill- 
conditioned matrices is of paramount importance for 
many problems, especially in the limit of large basis sets. 
Figure 5 illustrates the convergence for condition num- 
bers ranging from k = 10 1 to 10 4 . We use initial vectors 
based on Koopman's theorem in each case. The curves in 
Fig. 5 depict the convergence behavior for the lowest ex- 
citation energy for each condition number k, but similar 
patterns are observed for all of the first five eigenvalues 
in each case. The dashed (black) and dotted (red) lines 
show convergence patterns for the more well-conditioned 
systems, whereas the solid (green) and dash-dotted (blue) 
lines represent convergence for the less well-conditioned 
systems. The small blip observed at iteration 8 for the 



k = 10 3 system is an artifact of the error calculation 
because of a sign change relative to the reference value. 

The number of iterations required to reach the con- 
vergence, Err(n) < 10~ 12 , increases by almost an or- 
der of magnitude when the condition number is in- 
creased. This also indicates the potential improvement 
that could be reached by an efficient prcconditioner. For 
the two better-conditioned systems, the distribution of 
the smaller eigenvalues is more even, resulting in smooth 
curves and relatively rapid convergence. This pattern 
for the first eigenvalue is also observed in the more well- 
conditioned case in Figs. 3 and 4. 

In going from well- to ill-conditioned matrices, not only 
does the slope decrease, extending the number of itera- 
tions to convergence, but the morphology of the curves 
changes as well. A pronounced step/plateau pattern is 
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evident, particularly as the optimization proceeds. This 
behavior is typical of conjugate gradient schemes with 
clustered eigenvalues. 46 ' 84 ' 85 



C. Sensitivity to Numerical Error 

To probe the robustness of the molccular-orbital-frcc 
algorithm, we added noise of varying levels to a well- 
conditioned (ft = 10 2 ) system, as illustrated in Fig. 6. 
Again, we used a Koopman's theorem-based initial vector 
and observed the convergence behavior for the first eigen- 
value. Random noise in the range, r € ±[10~ 8 , 10 -4 ], 
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FIG. 6: Convergence of the first eigenvalue for k ~ 10 2 with 



random noise in the range r £ ±[10~ 
guesses based upon Koopman's theorem. 
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using initial 



provides a reasonable estimate of the induced errors 
we encounter in a typical linear scaling implementation, 
where small elements below some chosen numerical toler- 
ance are set to zero. We added this noise to every compo- 
nent of the newly formed vector (JLx) for each iteration 
of the inner conjugate gradient loop (the loop) in 
order to simulate the accumulation of numerical error as 
the calculation proceeds (see Fig. 1). 

We find that the algorithm is robust and stable with 
respect to numerical noise and that the error at con- 
vergence scales approximately linearly with the level of 
noise. The same behavior is also observed for more ill- 
conditioned systems, as shown in Fig. 7. 



DISCUSSION AND CONCLUSIONS 



Condition Number, k=10 
Initial Guess: Koopman's Theorem 
Convergence of Eigenvalue 1 



ffl 



60 
O 




T = 


io- 4 


t = 


io- 6 


" X = 


io- 8 


' ' T = 


IO" 10 




i I 



20 40 60 80 100 120 140 160 180 200 220 240 260 
Iteration (n) 



FIG. 7: Convergence of the first eigenvalue for k ~ 10 4 for 
random noise in the range, r £ ±[10 -8 ,10~ 4 ] using initial 
guess vectors derived from Koopman's theorem. 



ing) and numerical noise. This analysis clearly indi- 
cates a potential for reduced complexity calculations of 
large systems. However, there remain several open ques- 
tions: 1) The search space for the excitation spectrum 
corresponding to the dimensions of the Liouvillc oper- 
ator IL in the RPA eigenvalue problem scales quadrati- 
cally, 0(N 2 ), with system size. Unless the initial guess is 
highly accurate, we can expect the number of iterations 
required to reach convergence to increase with system 
size. This would obviate linear scaling complexity. 2) Un- 
fortunately, the construction of a highly accurate initial 
guess is computationally very expensive. For example, 
building the Koopman's guess would typically require a 
diagonalization of the Fockian which scales as 0(N 3 ). 
A reduced complexity technique for finding a good ini- 
tial guess or accurate preconditioning 52 ' 54 ' 56 ' 86,87 remains 
an unsolved, important challenge, though many efficient 
constructions may be possible. 3) The stability of the 
Wilkinson shift under sparse linear matrix algebra has 
not been fully investigated, though the stability under 
noisy conditions indicates that this is not a problem. 

In conclusion, the molecular-orbital-free scheme based 
upon a well-established variational characterization of 
the RPA excitation spectrum exhibits most of the neces- 
sary features required for an efficient linear scaling imple- 
mentation. While further work remains, we believe this 
technique will become highly valuable for determination 
of large-scale excited state properties. 



We have presented an algorithm for the variational 
characterization of the RPA eigenvalue spectrum, based 
on the Thouless functional and a non-linear conjugate 
gradient optimization in a molccular-orbital-free repre- 
sentation. We have analyzed the convergence with re- 
spect to initial guess, condition number (precondition- 
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